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LOW-THRUST TRAJECTORY OPTIMIZATION 
WITH SIMPLIFIED SQP ALGORITHM 


Nathan L. Parrish, and Daniel J. Scheerest 


The problem of low-thrust trajectory optimization in highly perturbed dynamics 
is a stressing case for many optimization tools. Highly nonlinear dynamics and 
continuous thrust are each, separately, non-trivial problems in the field of optimal 
control, and when combined, the problem is even more difficult. This paper de- 
scribes a fast, robust method to design a trajectory in the CRTBP (circular re- 
stricted three body problem), beginning with no or very little knowledge of the 
system. The approach is inspired by the SQP (sequential quadratic programming) 
algorithm, in which a general nonlinear programming problem is solved via a se- 
quence of quadratic problems. A few key simplifications make the algorithm pre- 
sented fast and robust to initial guess: a quadratic cost function, neglecting the 
line search step when the solution is known to be far away, judicious use of end- 
point constraints, and mesh refinement on multiple shooting with fixed-step inte- 
gration. 


INTRODUCTION 


Low-thrust, electric propulsion systems compare favorably to chemical propulsion systems be- 
cause the specific impulse (Isp) is generally an order of magnitude higher for electric propulsion. 
However, the use of low-thrust propulsion complicates mission design. While an impulsive trajec- 
tory can be fully represented with a finite, small number of variables, continuous-thrust trajectories 
are in principal of infinite dimensions. Thus, we must be clever to reduce the size of the continuous- 
thrust problem while preserving its accuracy. 


Most nonlinear optimization methods require an initial guess that is “close” to an optimal solu- 
tion, or at least “close” to a feasible solution. However, developing a good initial guess can be as 
difficult a problem as optimizing in the first place. In general, to find a good initial guess, we must 
make some assumption about the structure of the solution. Since there are many local optima to 
these problems, which one we arrive at is a function of the initial guess. The approach demonstrated 
in this paper is capable of finding optimal trajectories even when no initial guess is available, so it 
is not necessary to assume any particular structure for the solution. 


Another challenge that this paper overcomes 1s that general nonlinear programming (NLP) solv- 
ers can be very slow for specific problems. Some of the most widely-used NLP solvers include 
SNOPT!', IPOPT*, MATLAB fmincon, and KNITRO. Some examples of these NLP solvers in tra- 
jectory optimization include: NASA Goddard’s GMAT (General Mission Analysis Tool) and its 
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CSALT collocation tool; PSOPT (PseudoSpectral OPTimization) open-source collocation tool*>; 
NASA Goddard’s EMTG (Evolutionary Mission Trajectory Generator)°; NASA Johnson’s Coper- 
nicus Trajectory Design and Optimization System; and NASA JPL’s MALTO (Mission Analysis 
Low Thrust Optimizer)’. 


While these solvers are powerful and adept at solving a wide array of problems, the authors of 
the present work found that they can be very slow for low-thrust trajectory optimization — particu- 
larly as the problem size grows or as the linear assumptions internal to the optimizer break down. 
With the exception of IPOPT (which is open-source), high-performance NLP solvers are generally 
proprietary and each is set up as an “engineering black box”. When the algorithm encounters diffi- 
culties with an optimization problem, there is very little the user can do to improve the situation. 
The present research implements a simplified version of the SQP (sequential quadratic program- 
ming) algorithm in such a way that trajectory optimization problems are solved in a small fraction 
of the total number of iterations and thus, a small fraction of the total computation time. 


THE CIRCULAR RESTRICTED 3 BODY PROBLEM 


The circular restricted three body problem (CRTBP) model assumes that the spacecraft 1s mass- 
less in comparison to the two primaries (here, the Earth and Moon). The primaries orbit their bar- 
ycenter in a circular orbit. 


Dimensionless units are used so that all state variables are on the order of 1. One distance unit 
(DU) is equal to the mean distance between Earth and Moon, or 384747.962856037 km. The non- 
dimensional mass ratio is equal to 0.012150585609624. A synodic reference frame is used such 
that the Earth is fixed at the point [—, 0, 0]’, and the Moon is fixed at the point [1 — p, 0,0]. This 
reference frame is illustrated in Fig. 1. 
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Fig. 1. The Earth-Moon synodic reference frame, with libration points L; through Ls la- 
beled. Note that Earth and Moon are both fixed along the x-axis at (- 2) and (1 — p), respec- 
tively. 


The equations of motion for the system in this rotating frame are given by 
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The use of simplified dynamical models here is justified because the algorithm’s effectiveness 
would not be affected by higher-fidelity perturbations. 


THE OPTIMAL CONTROL PROBLEM 


We find it helpful to write the optimal control problem in different ways. Ultimately, the prob- 
lem we want to solve is: 


minimize: fpatn = J luce) II? dt (4) 


subject to: h(x, u,t)= 0 
2 (5) 
G(%U,t) <0 
When the exponent p in the objective function is equal to 1, we are finding the minimum fuel 
solution. When p = 2, we are finding what is commonly referred to as the “minimum energy” 
solution. The radius of convergence is much larger when p = 2, so we use that cost function ini- 
tially. In this paper, we will leave p = 2, but the present authors and others have shown previously 
how homotopy can be used to transition to the minimum fuel solution? ™. 


In this work, the equality constraints h(z, u, t) are used to enforce the problem dynamics and to 
constrain the endpoints. The inequality constraints g(%, u, t) are used to limit the thrust magnitude. 
We define the Lagrangian functional associated with the NLP as 


L(z,U, tA, fi) = f(#,U,t) +A-hG, t,t) + A: G3, 4, t). (6) 


Along with the KKT conditions constraining the Lagrange multipliers A and li, the NLP problem 
is now rewritten in the following form. 


minimize: L(%, i,t, A, ji) (7) 
subject to: h(x, u,t) = 0 
2770 CO m4 (8) 
g(%,u,t) < 0 


where £, h, and g are all nonlinear functions. 


SEQUENTIAL QUADRATIC PROGRAMMING 


Since it is, in general, hard or impossible to analytically solve a nonlinear problem, we need to 
convert the problem into some form that is possible to solve. The Sequential Quadratic Program- 
ming (SQP) algorithm is a robust option. The basic concept is to solve a series of quadratic sub- 


problems, each of which approximate the nonlinear problem at the current iteration’'*. The quad- 
ratic problem at each problem is constructed as follows: 


minimize: VL(Xp, Ux, C, iste) . OX + ~5X : HL (Xx, Ux, t, Desig) F OX (9) 


subject to: A(X, Ur, t) + Vh(X_, Ux, t): 8X = 0 ae 
G(Xq, Ug, t) + VG (Xp, liz, t) 5X <0 


where O;, are the variables about which the current iteration is linearized and 6X is the update to 
all of the optimization variables. After solving the quadratic problem at iteration k, the optimization 
variables are updated as 
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where @ is a Scalar in the range (0, 1] and is determined by performing a line search. Here, the line 
search is implemented by testing a range of a@ and choosing the value that minimizes the defect 
constraint violations. It was found that when solving with a poor initial guess, it 1s best to leave 
a = 1 for the first several iterations. Then, when the solution is “near” converging, perform the 
line search at every iteration. 


SECOND-ORDER CONIC PROGRAMMING 


Second-Order Conic Programming (SOCP) is a well-defined class of optimization problems for 
which many well-established, high-performance solvers exist. If the optimal control problem can 
be approximated as a SOCP, then we can solve the approximate problem quickly and efficiently. 
An SOCP problem consists of some set of the following: linear equality constraints; quadratic ine- 
quality constraints; linear objective; and quadratic objective. SOCP and QP (quadratic program- 
ming) are very similar, with the distinction that QP only allows inequality constraints to be up to 
linear, while SOCP allows inequality constraints to be up to quadratic. In the present implementa- 
tion, quadratic inequality constraints were used to limit the square of the thrust magnitude (which 
is identical to limiting the thrust magnitude). 


The JuMP open-source modeling language!’ was used in the Julia programming language!*. A 
great advantage of the JuMP modeling language is that it provides a common, user-friendly inter- 
face for a variety of solvers. Two solvers were used in this work: IPOPT (Interior Point OPTimizer), 
which is open-source, and Gurobi, which is free for academia. It was found that for these problems, 
Gurobi is some 20-40X faster, making it the clear choice when a license is available. JuMP can 
also call other SOCP solvers, but IPOPT and Gurobi were the only two that worked straight out of 
the box. 


MULTIPLE SHOOTING 


A trajectory is defined as a set of N nodes, where each node has position, velocity, mass, and 
thrust. For each iteration of the SQP algorithm, the relative time between each node is fixed, though 
the total time of flight is an optimization variable. For each pair of nodes, the state 1s propagated 
forward halfway from the first and backward halfway from the second. Continuous thrust is held 


constant for each node. The defect d in state at the midpoint is constrained to be zero. 
The Jacobian of the defect constraints with respect to all optimization variables is used to line- 


arize the dynamics constraints. Forward finite differencing is used to approximate each of the par- 
tial derivatives in the large matrix. The size of the Jacobian matrix grows with N*, where N is the 


number of nodes used. However, the number of nonzero elements grows with N. Since the Jacobian 
is Sparse, We Save computation time by computing only the partials that may be nonzero. 


The optimal control problem is cast in terms of an SOCP largely by enforcing a linearization of 
the dynamics constraints: 


k aR, (12) 


where d x 18 the vector of all defect constraints at the current iteration. Note that Eqn. (12) 1s equiv- 
alent to the more general Eqn. (10). The cost function of the original problem given in Eqn. (5) 1s 
truly quadratic, so that no approximation needs to be made. The objective can be left as-is. 


If the endpoints (initial and final nodes) and time of flight are held constant, then the dynamics 
constraints and objective function alone are enough to find optimal solutions. Adding these few 
extra variables to the optimization problem significantly increases its complexity, so a full section 
is devoted to it below. 


Fixed-step integration is used, for three reasons. First, when partial derivatives are computed 
via finite differencing, fixed-step integration yields much more consistent derivatives than adap- 
tive-step integration. Each adaptive-step integration adds a small amount of algorithmic noise be- 
cause the derivatives function is evaluated at slightly different times for a perturbed initial state. 
When this noise is divided by a small perturbation size in finite differencing, the noise is amplified 
and can in some cases be on the order of the derivative itself. Second, the integration 1s faster 
because the algorithm does not get hung up propagating near a gravity well. Finally, the amount of 
work done to propagate each node is identical, which is favorable for future implementation with 
parallel computing, either on multi-threaded CPU or on GPU. 


Mesh Refinement 


The reader may question: What about when the trajectory does go near a gravity well? Will the 
integration accuracy not degrade? To address this important concern, we implement mesh refine- 
ment, where nodes are added (or removed) such that the integration error remains within some 
fixed bounds. Runge-Kutta 78 integration is used where a 7" order polynomial is used to propagate 
the state, and the 8" order term is used to estimate the truncation error. Typically, an adaptive-step 
integrator would use the 8" order term to determine each integration step size. Here, we use a fixed 
step size and instead use the 8" order term to determine where more nodes are needed. 


We see in Figure 2 how a coarse initial solution can be refined into a more accurate solution. 
The mesh refinement operation is fast, requiring on the order of tens of milliseconds. In addition to 
improving the numerical integration accuracy, mesh refinement tends to assist with the lineariza- 
tion of the problem. Portions of a trajectory which are deep in a gravity well have quickly-changing 
dynamics, and changes to those nodes have a great effect on the rest of the trajectory. Mesh refine- 
ment addresses the accuracy and sensitivity of these nodes simultaneously. 
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Figure 2. An exaggerated example of the mesh refinement operation. Here, a transfer from 
a DRO to an Lz halo orbit was first solved with just 10 nodes in (a), then refined to (b). The 
plots are in the Earth-Moon rotating frame. Control is held constant for each node. 


ENDPOINT CONSTRAINTS 


In the simplest terms, we are solving a two-point boundary value problem. It was found that 
using fully-constrained endpoints (i.e. fixed position, velocity, and time) made the problem much 
easier to solve because the linear approximation breaks down for the endpoints sooner than for the 
other optimization variables. However, for many mission design scenarios, we need to be able to 
optimize the endpoints and time of flight in addition to the intermediary states and controls. When 
the end states are instead constrained to lie anywhere on an orbit, the 6-state constraint 1s reduced 
to a 5-state constraint. If the initial and final orbits are Keplerian, then it is natural to describe the 
state as a Set of orbital elements where, for instance, true anomaly is allowed to drift freely. If the 
end states are instead constrained to lie on periodic three-body orbits, then there is no equivalent 
set of only 5 elements to constrain. Instead, we must use a linear or quadratic Taylor series expan- 
sion of the three-body orbit with respect to the “true anomaly”. We will use the angle T to be similar 
to true anomaly and represent the angular distance traveled around a general periodic orbit. 


For any arbitrary dynamics and endpoint, we define the function g(t) to be the actual endpoint 
state in Cartesian position and velocity as a function of the angle tT. If g represents a periodic orbit, 
then g(0) = q(27). As implemented in this work, g(t) is defined by a set of states given at 100 
equally-spaced values of T, read from a text file. The endpoints are interpolated from the set of 100 
states given. 


The Taylor series expansion of the endpoints are found as follows. First, compute the partial 
derivatives of the endpoint with respect to Tt via finite differencing. 
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Note that with only three evaluations of g(t), we can estimate both the first and second derivatives. 
Then, we can form an approximate, quadratic form of the endpoint orbit from the first two terms 
of the Taylor series. 
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In SOCP, we can use up to linear equality constraints and up to a quadratic objective, but we 
cannot use a quadratic equality constraint. It was found that using only a linear equality constraint, 
the optimization algorithm would tend to bounce around the optimal T, never quite converging. We 
developed two solutions to this problem. The first solution is to introduce a quadratic term to the 
objective function which measures the distance between the quadratic expansion and the linear 
expansion. Adding this quadratic cost to the path cost yields the following objective function: 
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where f is a scaling term. If 6 is too small, then the solution bounces around the optimal T indefi- 
nitely. If 6 is too large, then the solution will converge quickly (but prematurely) on a sub-optimal 
value of T. 


The other approach to help the endpoints converge is simpler: only allow T to vary on some 
iterations. For instance, constrain T to be fixed on the even iterations and allowed to vary on the 
odd iterations. When the solution is “nearly” converged, then we can safely fix the endpoints and 
know that we are at a “nearly” optimal solution. As mentioned before, once the endpoints are fixed, 
the problem tends to converge the rest of the way within a few iterations. Both approaches for 
optimizing the endpoints are somewhat ad-hoc, and admittedly are not guaranteed to converge to 
the optimal endpoints. However, it was found that the endpoints do consistently converge to similar 
values, giving us confidence that the algorithm is effective. It is important to note that for low- 
thrust propulsion, the endpoint orbits are departed or arrived at gradually, making the exact value 
of the t variables less important than for impulsive-burn trajectories. For instance, two solutions 
could be identical if t and time of flight are adjusted together — one solution would simply have a 
coasting period at the end. 
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Figure 3. The linear and quadratic approximations of an Earth-Moon Lz halo orbit. Both 
the linearized endpoint and the quadratic endpoint are tangent to the halo orbit at one point, 
but the quadratic endpoint remains close to the halo orbit for a longer time. 


Constraining the first and last nodes of the trajectory to lie on a linear expansion of a 3-body 
orbit was found to be effective. The endpoints are constrained to lie anywhere along a 6-dimen- 
sional line (3 dimensions for position, 3 for velocity). Figure 3 shows an example of a linear and 
quadratic expansion of an L2 halo orbit. 
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Figure 4. Linear expansions of the modified equinoctial elements representation of a halo 
orbit. The black dotted line is the linear expansion. The approximation holds well near posi- 
tion (a), less well at position (b), and breaks down at position (c) due to the singularity at 180° 
inclination. 


As mentioned earlier, two-body orbits can easily be parameterized with orbital elements, mak- 
ing it natural to constrain the 5 constant elements and optimize the final element. We explored 
using Modified Equinoctial Elements relative to the Moon to describe a halo orbit and found that 
the orbital element set breaks down due to singularities. Hintz compiled a list of 22 orbital element 
sets!> and found that each one has at least one singularity. In the case of the Modified Equinoctial 
Elements, the singularity is at inclination of 180°. For two-body orbits, such a singularity is easy to 
avoid, but even a simple three-body orbit can cross every possible singularity. Figure 4 shows how, 
in some cases, a linear expansion of a halo orbit represented with Modified Equinoctial Elements 
is more accurate than a linear expansion of the same orbit represented with Cartesian position and 
velocity. In other cases, however, the orbital elements are far less accurate because of the singular- 


ity. 


EXAMPLES & RESULTS 


The greatest benefits from the approach described in this paper are the robustness to poor initial 
guess and the small number of iterations required to converge on solutions. To demonstrate these 
qualities, we provide the worst-possible initial guess: random noise. The states guesses are pulled 
from a random number generator, and the controls guess is zero. Even with such a poor initial 
guess, the algorithm will still converge quickly. 
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Figure 5. Progression from random noise in (a) to a converged solution in (d) with 19 iter- 
ations, for a DRO-to-DRO transfer with fixed endpoints and time of flight. The first DRO (in 
blue) has an x-crossing at 0.9 DU. The second DRO (in orange) has an x-crossing at 0.8 DU. 
Earth and Moon are plotted to scale in the synodic reference frame. The canonical distance 
unit (DU) is used where 1 DU is equal to the average distance between Earth and Moon. 


This ability to rapidly converge on an optimal solution is demonstrated in Fig. 5. In almost any 
situation, it is possible to come up with a better initial guess than random noise. Thus, this result 
should be considered a worst case. If we really have no knowledge of what the solution should look 
like, even a random guess is close enough to converge. It should be noted that this approach works 
well for transfers up to about two complete revolutions. Beyond that, the linearity assumptions in 
the dynamics constraints break down, and the line search step becomes increasingly necessary. 
Future work will seek to improve how the algorithm handles multiple-revolution transfers. 


We will now show some more example solutions for transfers between DRO’s and halo orbits. 
All of these solutions began with a random initial guess. Thrust was unconstrained in most cases, 
but the objective function (integral of the control squared) encouraged low-thrust solutions. From 
prior work, we know that each of these can be converted into a nearby minimum fuel solution with 
thrust magnitude requirements less than or equal to the maximum thrust stated. 
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Figure 6. An example transfer from one L2 halo orbit to another, in the Earth-Moon sys- 
tem. The first halo orbit is in blue and the second in orange. This transfer requires 26 days 
and an acceleration of 3.5E-5 m/s” (equivalently, 35 mN thrust for a 1000 kg spacecraft). 
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Figure 7. Example transfer from a DRO to another DRO. The first DRO is in blue and 
crosses the x-axis at 0.8; the second DRO is on orange and cross the x-axis at 0.9. This transfer 
requires 15 days and an acceleration of 1.7E-4 m/s” (equivalently, 170 mN for a 1000 kg space- 
craft). 
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Figure 8. An example transfer from a DRO to an Lz halo orbit. The DRO is in blue and 
the halo orbit in orange. This transfer requires 29 days and an acceleration of 2.8E-4 m/s” 
(equivalently, 280 mN for a 1000 kg spacecraft). 
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Figure 9. An example transfer from an Lz halo orbit to a DRO (the reverse of Fig. 7). The 
halo orbit is in blue and the DRO in orange. This transfer requires 28 days. The max control 
acceleration was constrained to be 2.8E-4 m/s” to match the prior example. 


Figures 6-9 provide some example optimal solutions. Note that the thrust levels required are 
within the range of what current electric propulsion technology can provide. Future work will ex- 
plore many more optimal solutions between more 3-body orbit types. 


Computation Time 


Based on the authors’ experience with other tools capable of solving similar problems, it was 
found that the present implementation compares favorably. Future work will include a detailed 
comparison of this tool with others. For now, we have results for the typical time required by the 
present research to solve the example problems shown. 
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Table 1. Comparison of computation time for each part of an iteration 


Iteration segment Typical time 
Set up SOCP problem | 0.2 —0.5 seconds 
Solve SOCP problem | 0.2 — 0.5 seconds 


Line search 0.2 — 0.5 seconds 


As Table | shows, the computational load is well balanced across the three main segments of 
an iteration: setting up the SOCP problem, solving the SOCP problem, and performing the line 
search. Each requires about 0.2 seconds for 40 nodes or about 0.5 seconds for 100 nodes. The SOCP 
setup time largely consists of computing the Jacobian of dynamics constraints with respect to states 
and controls. The Jacobian calculation can potentially be accelerated by computing the Jacobian in 
parallel on the GPU. Solving the SOCP problem is done by the Gurobi solver, which is already 
highly optimized and one of the fastest solvers available. We do not expect that any improvement 
can be made to solving the SOCP problem, and such work is considered outside the scope of the 
present research. The line search can be accelerated by using a more efficient method to find the 
best a. By accelerating the SOCP problem setup and the line search, we expect that the total process 
can be accelerated by a factor of two still. 


The exact number of iterations required varies depending on the quality of the initial guess. The 
number if iterations stated is a typical value from running these transfers several times with differ- 
ent random initial guesses. The total number of iterations required for a problem with fixed end- 
points is about 10-20, meaning the total solution time is roughly 2-10 seconds starting from a ran- 
dom initial guess. With a close initial guess, the number of iterations required is typically 2-5, 
meaning 1-2 seconds to a converged, optimal solution. When the endpoints are optimized, the al- 
gorithm typically requires about twice as many iterations. 


CONCLUSION 


In comparison to the traditional approach of plugging the problem into a “black-box” NLP 
solver, the methods shown converge even when given no knowledge of the solution at all. It was 
found that the only piece of information that the user needs to provide is a rough guess for the time 
of flight, as the transfer time guess will dictate which set of local solutions the algorithm could 
converge on. This robustness to initial guess is a compelling feature, as three-body orbit transfers 
are challenging to design with intuition alone. Of course, if a high-quality initial guess is available, 
the methods shown are still valid. 


We have shown that endpoints can be efficiently constrained to lie on 3-body repeating orbits, 
and that time of flight can be optimized as well. When optimizing the endpoints, we must make a 
trade between converging quickly on sub-optimal endpoints or converging more slowly on end- 
points that are arbitrarily close to optimal. It 1s easy for the mission design engineer to adjust this 
trade based on the problem at hand. 


The biggest limitation to the algorithm at this point is that multi-revolution transfers (greater 
than 2 revolutions) do not work nearly as well. This restriction comes in because the relationship 
between node | and node N becomes increasingly nonlinear as the angular distance grows. Trans- 
fers with more than about 1.5 complete revolutions generally require the line search to improve 
convergence. 
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Future work includes: Comparison of this algorithm with other established tools; improvements 
to how multiple-revolution transfers are handled; parallelization of the Jacobian computation; i1n- 
creased efficiency for the line search; and optimization of many more trajectories between a variety 
of 3-body orbits. 
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